Fast stable REML and ML estimation of semiparametric GLMs

Recent work by Reiss and Ogden (2009) provides a theoretical basis for sometimes preferring restricted maximum likelihood (REML) to generalized cross validation (GCV) for smoothing parameter selection in semiparametric regression. However, existing REML or marginal likelihood (ML) based methods for semiparametric GLMs use iterative REML/ML estimation of the smoothing parameters of working linear approximations to the GLM. Such indirect schemes need not converge, and fail to do so in a non-negligible proportion of practical analyses. By contrast, very reliable prediction error criteria smoothing parameter selection methods are available, based on direct optimization of GCV, or related criteria, for the GLM itself. Since such methods directly optimize properly defined functions of the smoothing parameters, they have much more reliable convergence properties. This article develops the first such method for REML or ML estimation of smoothing parameters. A Laplace approximation is used to obtain an approximate REML or ML for any GLM, which is suitable for efficient direct optimization. This REML/ML criterion requires that Newton-Raphson, rather then Fisher scoring, be used for GLM fitting, and a computationally stable approach to this is proposed. The REML or ML criterion itself is optimized by a Newton method, with the required derivatives obtained by a mixture of implicit differentiation and direct methods. The method will cope with numerical rank deficiency in the fitted model, and in fact provides a slight improvement in numerical robustness on the method of Wood (2008) for prediction error criteria based smoothness selection. Simulation results suggest that the new REML and ML methods offers some improvement in mean square error performance relative to GCV/AIC in most cases, without the small number of severe undersmoothing failures to which AIC and GCV are prone. This is achieved at the same computational cost as GCV/AIC. The new approach also eliminates the convergence failures of previous REML/ML based approaches for penalized GLMs, and usually has lower computational cost than these alternatives. Example applications are presented in adaptive smoothing, scalar on function regression and generalized additive model (GAM) selection.